Introduction

To facilitate transparency, the following document shows code behind the results of the manuscript. The display items are ordered according to their appearance in the main text. Note that Table R just summarises the model outputs used to extract the repeatability estimates. Additional figures, used in the responses to the reviewers are also shown.

The code is displayed upon clicking the “show” icon at the top right, above each display item.

When using this content PLEASE CITE the preprint and this repository: Bulla & Sankar, 2025. Supporting information for ‘No support for honest signalling of male quality in zebra finch song’. GitHub https://github.com/MartinBulla/rebuttal_alam_2024.

Questions about the outpus can be directed to and about the UMAP generation to

Repository: files & folders

Supportin information, including code: the current html document with all display items and code.

Data:

  • Fig_4c.csv: source data provided by Alam et al. (2024, Nature); columns: tutor - unique tutor ID, pupil - unique pupil ID for the given tutor, tutor_path - song path length of the tutor, pupil_path - song path length of the pupil
  • Dat_path_length.csv: 20 UMAP iterrations generated from 31 tutored birds of Alam et al.’s Fig 2c 2, using Alam et al.’s methods to:
  • Fig_2b.csv: source data provided by Alam et al. (2024, Nature) for their Fig. 2b; columns: bird - unique bird ID, 1-15 - columns names 1-15 indicate 15 iterations and contain path-length estimates for each
  • rebuttal_fig_2.csv: source data provided by Alam et al. (2024, Nature) for their Fig. 3c and data measured by us from their Extended Data Fig. 6 (green versus black bars during playback); columns: trial_id - unique trial ID, side_bias - was the long-path song during the trial played in the preferred arm or non-preferred arm by the female during the baseline period?, song_pair - unique song-pair ID; long - path length of the long path song, short - path length of the short path song, pre - the percentage of pre-trial time spent in the arm where long path song would be played during the trial, trial - the percentage of trial time spent in the arm with long path song, post - the percentage of post-trial time spent in the arm where long path song was played during the trial, trial_long - trial time spent in the arm with long path song, trial_short - trial time spent in the arm with short path; the dataset was used during the initial submission
  • Fig_3c_ED_Fig_6.csv: source data provided by Alam et al. (2024, Nature) for their Fig. 3c and data measured by us from their Extended Data Fig. 6 (green versus black bars during playback); columns: trial_id - unique trial ID, side_bias - was the long-path song during the trial played in the preferred arm or non-preferred arm by the female during the baseline period?, song_pair - unique song-pair ID; long - path length of the long path song, short - path length of the short path song, pre - the percentage of pre-trial time spent in the arm where long path song would be played during the trial based on source data from Fig. 3c, trial - the percentage of trial time spent in the arm with long path song based on source data from Fig. 3c, post - the percentage of post-trial time spent in the arm where long path song was played during the trial based on source data from Fig. 3c, pre_ed - the percentage of pre-trial time spent in the arm where long path song would be played during the trial as measured from Extended Data Fig. 6, trial_ed - the percentage of trial time spent in the arm with long path song as measured from Extended Data Fig. 6, post_ed - the percentage of post-trial time spent in the arm where long path song was played during the trial as measured from Extended Data Fig. 6, comments note the discrepencies between source data of Fig. 3c and Extended Data Fig. 6
  • match_labels_output.csv: contains matched bird ids and names from the Fig_4c.csv and Fig. 2b
  • Dat_path_length.csv generated by the procedure described here and used to add syllable counts to Fig_2b.csv data; columns: f2b_bird_id - unique bird ID from Fig. 2b, f2b_bird_name - unique bird name from Fig. 2b, superset_bird_id - unique bird ID from Fig. 4c, Alam et al’s file tut.zip, superset_bird_name - unique bird name from Fig. 4c, Alam et al’s file tut.zip, superset_bird_cleaned_name - Fig. 4c bird names cleaned to match the Fig. 2b names, superset_bird_type - was the bird isolated or tutored, superset_bird_nsyll - number of unique syllables in a given song
  • Dat_rnorm.Rdata: random data from the Fig 1 generated in the Fig_point_1.r script

R-scripts used in the analysis,❗run relative to the project’s root directory, with the same folder structure as the repository:

For the first submission:

For the revision

Output: contains the generated figures in PNG format

LICENSE: terms of reuse

Code to load tools & data
# load R-packages
require(arm)
require(cowplot)  
require(data.table)
require(foreach)     
require(ggplot2)
require(ggpmisc)
require(ggpubr)
require(grafify)
require(grid)
require(gridExtra)
require(kableExtra)
require(legendry) 
require(pander) 
require(patchwork)
require(rptR)
require(tidyr)

# constants
set.seed(5) # for simulations
nsim = 5000
g_r_x = 35
orig_ = '#D43F3AFF'
furt_ = '#46B8DAFF'#"#357EBDFF"#
col_p1 = 'darkgrey'#'darkgrey'
col_p = 'grey60'#'darkgrey'
col_pc = 'grey10'
point_out = 'grey30'
point_fill = "#46B8DAFF"
col_R = "#D43F3AFF"
col_l = '#D43F3AFF'
col_ <- c("#46B8DAFF", "#EEA236FF")
#c("#357EBDFF", "#D43F3AFF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#9632B8FF", "#9632B8FF")[7:1]
#col_2 <- c("#357EBDFF", "#D43F3AFF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#9632B8FF", "#9632B8FF") # require(scales);show_col(col_2)
#col_3 <- c("#46B8DAFF","#5CB85CFF", "#EEA236FF") # require(scales);show_col(col_2)
per_ = 0.05
scale_size = 0.352778
inch = 0.393701
force_0_yes = FALSE

# function for model outpus
m_out = function(model = m, name = "define", 
        type = "mixed", 
        dep = "define", fam_ = 'Gaussian',
        round_ = 2, nsim = 5000, aic = FALSE, save_sim = here::here('Data/model_sim/'), back_tran = FALSE, perc_ = 1, R2 = FALSE){
          # perc_ 1 = proportion or 100%
        bsim = sim(model, n.sim=nsim)  
        
        if(save_sim!=FALSE){save(bsim, file = paste0(save_sim, name,'.RData'))}
       
        if(type != "mixed"){
          v = apply(bsim@coef, 2, quantile, prob=c(0.5))
          ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975)) 

          if(back_tran == TRUE & fam_ == "binomial"){
           v = perc_*plogis(v)
           ci = perc_*plogis(ci)
           }
          if(back_tran == TRUE & fam_ == "binomial_logExp"){
                v = perc_*(1-plogis(v))
                ci = perc_*(1-plogis(ci))
                ci = rbind(ci[2,],ci[1,])
               }

          if(back_tran == TRUE & fam_ == "Poisson"){
           v = exp(v)
           ci = exp(ci)
          }

         oi=data.frame(type='fixed',effect=rownames(coef(summary(model))),estimate=v, lwr=ci[1,], upr=ci[2,])
          rownames(oi) = NULL
          oi$estimate_r=round(oi$estimate,round_)
          oi$lwr_r=round(oi$lwr,round_)
          oi$upr_r=round(oi$upr,round_)
          if(perc_ == 100){
           oi$estimate_r = paste0(oi$estimate_r,"%")
           oi$lwr_r = paste0(oi$lwr_r,"%")
           oi$upr_r = paste0(oi$upr_r,"%")
          }
         x=data.table(oi[c('type',"effect", "estimate_r","lwr_r",'upr_r')]) 
       
        }else{
         v = apply(bsim@fixef, 2, quantile, prob=c(0.5))
         ci = apply(bsim@fixef, 2, quantile, prob=c(0.025,0.975)) 

         if(back_tran == TRUE & fam_ == "binomial"){
          v = perc_*plogis(v)
          ci = perc_*plogis(ci)
         }
          if(back_tran == TRUE & fam_ == "binomial_logExp"){
                v = perc_*(1-plogis(v))
                ci = perc_*(1-plogis(ci))
                ci = rbind(ci[2,],ci[1,])
               }

          if(back_tran == TRUE & fam_ == "Poisson"){
            v = exp(v)
            ci = exp(ci)
         }

        oi=data.table(type='fixed',effect=rownames(coef(summary(model))),estimate=v, lwr=ci[1,], upr=ci[2,])
            rownames(oi) = NULL
            oi[,estimate_r := round(estimate,round_)]
            oi[,lwr_r := round(lwr,round_)]
            oi[,upr_r :=round(upr,round_)]
            if(perc_ == 100){
             oi[,estimate_r := paste0(estimate_r,"%")]
             oi[,lwr_r := paste0(lwr_r,"%")]
             oi[,upr_r := paste0(upr_r,"%")]
            }
         oii=oi[,c('type',"effect", "estimate_r","lwr_r",'upr_r')] 
        
         l=data.frame(summary(model)$varcor)
         l = l[is.na(l$var2),]
         l$var1 = ifelse(is.na(l$var1),"",l$var1)
         l$pred = paste(l$grp,l$var1)

         q050={}
         q025={}
         q975={}
         pred={}
         
         # variance of random effects
         for (ran in names(bsim@ranef)) {
           ran_type = l$var1[l$grp == ran]
           for(i in ran_type){
            q050=c(q050,quantile(apply(bsim@ranef[[ran]][,,ran_type], 1, var), prob=c(0.5)))
            q025=c(q025,quantile(apply(bsim@ranef[[ran]][,,ran_type], 1, var), prob=c(0.025)))
            q975=c(q975,quantile(apply(bsim@ranef[[ran]][,,ran_type], 1, var), prob=c(0.975)))
            pred= c(pred,paste(ran, i))
            }
           }
         # residual variance
         q050=c(q050,quantile(bsim@sigma^2, prob=c(0.5)))
         q025=c(q025,quantile(bsim@sigma^2, prob=c(0.025)))
         q975=c(q975,quantile(bsim@sigma^2, prob=c(0.975)))
         pred= c(pred,'Residual')

         ri=data.table(type='random',effect=pred, estimate_r=round(100*q050/sum(q050)), lwr_r=round(100*q025/sum(q025)), upr_r=round(100*q975/sum(q975)))
           
         ri[lwr_r>upr_r, lwr_rt := upr_r]
         ri[lwr_r>upr_r, upr_rt := lwr_r]
         ri[!is.na(lwr_rt), lwr_r := lwr_rt]
         ri[!is.na(upr_rt), upr_r := upr_rt]
         ri$lwr_rt = ri$upr_rt = NULL

         ri[,estimate_r := paste0(estimate_r,'%')]
         ri[,lwr_r := paste0(lwr_r,'%')]
         ri[,upr_r := paste0(upr_r,'%')]
        
        x = data.table(rbind(oii,ri))
        }
        
        x[1, model := name]                                                                
        x[1, response := dep]                                                                
        x[1, error_structure := fam_]      
        N = length(resid(model))                                                          
        x[1, N := N ]                                                                

        x=x[ , c('model', 'response', 'error_structure', 'N', 'type',"effect", "estimate_r","lwr_r",'upr_r')] 

        if (aic == TRUE){   
            x[1, AIC := AIC(update(model,REML = FALSE))] 
            }
        if (aic == "AICc"){
            aicc = AICc(model)
            x[1, AICc := aicc] 
        }
        if(type == "mixed" & nrow(x[type=='random' & estimate_r =='0%'])==0 & R2 == TRUE){
          R2_mar = as.numeric(r2_nakagawa(model)$R2_marginal)
          R2_con = as.numeric(r2_nakagawa(model)$R2_conditional)
          x[1, R2_mar := R2_mar]
          x[1, R2_con := R2_con]
         }
        x[is.na(x)] = ""
        return(x)
      }

Point 1 - statistical artifact

# Fig 1, top panels - random data
h = data.table(
x = rnorm(1000, mean = 35, sd = 4),
y = rnorm(1000, mean = 35, sd = 4)
)
#save(file='Data/Dat_rnorm.Rdata', h)
load(file=here::here('Data/Dat_rnorm.Rdata')) # loads the same sample as in the plot

g1=
ggplot(h, aes(x = x, y = y)) + 
    geom_point(col = col_p1, size = 1.5, alpha = 0.6) +
    stat_poly_line(se = FALSE, col = col_l) +
    stat_cor(cor.coef.name = "r", aes(label = after_stat(r.label)),  col = col_R, r.accuracy = 0.1, label.x = g_r_x, label.y = 52-(52-20)*per_ , hjust = 0.5, cex = 3) + 
    coord_cartesian(xlim = c(20, 50), ylim = c(20, 52)) +
    scale_x_continuous(breaks = seq(20,50, by=10), labels = seq(20,50, by=10), expand = c(0,0))+
    scale_y_continuous(breaks = seq(20,50, by=10), labels = seq(20,50, by=10), expand = c(0,0)) +
    theme_bw() + 
    theme(axis.ticks = element_blank(),
        panel.grid.major = element_blank(),  # Remove major grid lines
        panel.grid.minor = element_blank() )
    #stat_poly_eq() 
g2=  
ggplot(h, aes(x = x, y = y-x)) + 
    geom_point(col = col_p1, size = 1.5, alpha = 0.6) +#geom_point(fill = col_p, pch = 21, size = 1.5) +#geom_point(col = col_p) +
    stat_poly_line(se = FALSE, col = col_l) +
    stat_cor(cor.coef.name = "r", aes(label = after_stat(r.label)),  col = col_R, r.accuracy = 0.1, label.x = g_r_x, label.y = 20-36*per_, hjust = 0.5, cex = 3) +
    coord_cartesian(xlim = c(20, 50), ylim = c(-20, 20)) +
    scale_x_continuous(breaks = seq(20,50, by=10), labels = seq(20,50, by=10), expand = c(0,0))+
    scale_y_continuous(breaks = seq(-20,20, by=10), labels = seq(-20,20, by=10),expand = c(0,0)) +
    theme_bw() + 
    theme(axis.ticks = element_blank(),
        panel.grid.major = element_blank(),  # Remove major grid lines
        panel.grid.minor = element_blank() )

g3=  
ggplot(h, aes(x = x, y = y+x)) + 
    geom_point(col = col_p1, size = 1.5, alpha = 0.6) +#geom_point(fill = col_p, pch = 21, size = 1.5) +
    stat_poly_line(se = FALSE, col = col_l) +
    stat_cor(cor.coef.name = "r", aes(label = after_stat(r.label)),  col = col_R, r.accuracy = 0.1, label.x = g_r_x, label.y = 90-(90-50)*per_, hjust = 0.5, cex = 3) + 
    coord_cartesian(xlim = c(20, 50), ylim = c(50, 90)) +
    scale_x_continuous(breaks = seq(20,50, by=10), labels = seq(20,50, by=10), expand = c(0,0))+
    scale_y_continuous(breaks = seq(50,90, by=10), labels =seq(55,95, by=10),expand = c(0,0)) +
    theme_bw() + 
    theme(axis.ticks = element_blank(),
        panel.grid.major = element_blank(),  # Remove major grid lines
        panel.grid.minor = element_blank() )

    #g1|g2|g3

    #ggsave(here::here('Output/Fig_1.png'), g1|g2|g3, width = 12*2, height = 4*2, units = 'cm')

# Fig 1, bottom panels - Alam et al's data from Fig. 4c
hh = fread(here::here('Data/Fig_4c.csv'))
hh[, pupil_path_full:=tutor_path+pupil_path]
hh[, tutor_minus_pupil:=pupil_path]

gg1=
ggplot(hh, aes(x =tutor_path, y = pupil_path_full)) + 
    geom_point(col = col_p, size = 2.2, alpha = 0.8) +
    scale_color_grafify(palette = 'fishy')+#, ColSeq = FALSE)+
    #annotate("text", x=g_r_x, y=50-(50-20)*per_*4, label= "Color indicates a tutor", cex = 2.75, hjust = 0.5)+ 
    stat_poly_line(se = FALSE, col = col_l) +
    stat_cor(cor.coef.name = "r", aes(label = after_stat(r.label)),  col = col_R, r.accuracy = 0.1, label.x = g_r_x, label.y = 50-(50-20)*per_, hjust = 0.5, cex = 3)+ 
    coord_cartesian(xlim = c(20, 50), ylim = c(20, 50)) +
    xlab('Tutor path length') + ylab('Pupil path length') +
    scale_x_continuous(breaks = seq(20,50, by=10), labels = seq(20,50, by=10), expand = c(0,0))+
    scale_y_continuous(breaks = seq(20,50, by=10), labels = seq(20,50, by=10), expand = c(0,0))+
    theme_bw() + 
    theme(  legend.position="none",
            axis.ticks = element_blank(),
        panel.grid.major = element_blank(),  # Remove major grid lines
        panel.grid.minor = element_blank() )

gg2=  
ggplot(hh, aes(x =tutor_path, y = pupil_path_full-tutor_path)) + 
    geom_point(col = col_p, size = 2.2, alpha = 0.8) +
    scale_color_grafify(palette = 'fishy')+#, ColSeq = FALSE)+
    stat_poly_line(se = FALSE, col = col_l) +
    stat_cor(cor.coef.name = "r", aes(label = after_stat(r.label)),  col = col_R, r.accuracy = 0.1, label.x = g_r_x, label.y = 20-(20+20)*per_, hjust = 0.5, cex = 3)+ 
    coord_cartesian(xlim = c(20, 50), ylim = c(-20, 20)) +
    xlab('Tutor path length') + ylab('Pupil minus tutor path length') +
    scale_x_continuous(breaks = seq(20,50, by=10), labels = seq(20,50, by=10), expand = c(0,0))+
    scale_y_continuous(breaks = seq(-20,20, by=10), labels =seq(-20,20, by=10),expand = c(0,0))+
    theme_bw() + 
    theme(  legend.position="none",
            axis.ticks = element_blank(),
        panel.grid.major = element_blank(),  # Remove major grid lines
        panel.grid.minor = element_blank() )
 
gg3=  
ggplot(hh, aes(x =tutor_path, y = tutor_path+pupil_path_full)) + 
    geom_point(col = col_p, size = 2.2, alpha = 0.8) +
    scale_color_grafify(palette = 'fishy')+#, ColSeq = FALSE)+
    stat_poly_line(se = FALSE, col = col_l) +
    stat_cor(cor.coef.name = "r", aes(label = after_stat(r.label)),  col = col_R, r.accuracy = 0.1, label.x = g_r_x, label.y = 90-(90-50)*per_, hjust = 0.5, cex = 3)+ 
    coord_cartesian(xlim = c(20, 50), ylim = c(50, 90)) +
    xlab('Tutor path length') + ylab('Pupil + tutor path length')+
    scale_x_continuous(breaks = seq(20,50, by=10), labels = seq(20,50, by=10), expand = c(0,0))+
    scale_y_continuous(breaks = seq(50,90, by=10), labels =seq(50,90, by=10),expand = c(0,0))+
    theme_bw() + 
    theme(  legend.position="none",
            axis.ticks = element_blank(),
        panel.grid.major = element_blank(),  # Remove major grid lines
        panel.grid.minor = element_blank()  )

    #gg1|gg2|gg3

    #ggsave(here::here('Output/Fig_2_v2.png'), gg1|gg2|gg3, width = 12*2, height = 4*2, units = 'cm')

# combine
combine_1 =  (g1|g2|g3) & plot_annotation(subtitle = "Random data; n = 1,000") #& theme(plot.title = element_text(hjust = .5))
combine_2 =  (gg1|gg2|gg3) & plot_annotation(subtitle = "Alam et al. 2024 data; n = 17 pupils and 8 tutors ") #& theme(plot.title = element_text(hjust = .5))

ggsave(here::here('Output/rev_Fig_1_width-137mm.png'), wrap_elements(combine_1) / wrap_elements(combine_2), width = 10.5*2, height = 4*2*2, units = 'cm')

wrap_elements(combine_1) / wrap_elements(combine_2)

Figure 1 | Illustration of illusory relationships when an x-variable is included in both axes. The top panels are based on the 1,000 randomly sampled values of x and y, the bottom panels depict data from Fig. 4c of Alam et al. 2024. In both the top and bottom panels, dots represent individual data points, with colour in the bottom panels indicating individual tutors. Lines represent ordinary least-square regressions. r denotes a Pearson’s correlation coefficient. The left panels highlight the absence of relationships in the data. The middle panels show negative relationships, and the right panels show positive relationships, both arising from including an x-variable also in the y-variable.

Output used in the replies to the reviewer comments

# replies to the referees - residuals from the line of idenity
## point A
g0 = 
ggplot(h, aes(x = x, y = y)) + 
    geom_segment(aes(xend = x, yend = x), color = "skyblue") +    # Residual lines
    geom_point(col = col_p) + # raw data
    geom_abline(intercept = 0, slope = 1, color = col_l) +  # Line of identity (y = x)
    labs(title = "Residuals from Line of Identity", subtitle = "Random data; n = 1000") +  # Labels
    coord_cartesian(xlim = c(20, 50), ylim = c(20, 52)) +
    scale_x_continuous(breaks = seq(20,50, by=10), labels = seq(20,50, by=10), expand = c(0,0))+
    scale_y_continuous(breaks = seq(20,50, by=10), labels = seq(20,50, by=10), expand = c(0,0))+
    theme_minimal()   

gg0 = 
ggplot(hh, aes(x = tutor_path, y = pupil_path_full)) + 
    geom_segment(aes(xend = tutor_path, yend = tutor_path), color = "skyblue") +    # Residual lines
    geom_point(col = col_p) + # raw data
    geom_abline(intercept = 0, slope = 1, color = col_l) +  # Line of identity (y = x)
    labs(title = "", subtitle = "Alam et al. 2024 data; n = 17 pupils and 8 tutors", xlab = "Tutor path length", ylab = "Pupil path length") +  # Labels
    coord_cartesian(xlim = c(20, 50), ylim = c(20, 52)) +
    scale_x_continuous(breaks = seq(20,50, by=10), labels = seq(20,50, by=10), expand = c(0,0))+
    scale_y_continuous(breaks = seq(20,50, by=10), labels = seq(20,50, by=10), expand = c(0,0))+
    theme_minimal() 

#ggsave(here::here('Output/rev_Reply_1.png'), grid.arrange(g0, gg0, ncol = 2), width = 14.7, height = 8, units = 'cm')
grid.arrange(g0, gg0, ncol = 2)

Responses Figure 1 | Illustration of the residual from the line of identity. Grey point represent individual observations (from Fig. 1), red line the line of identity betwee x and y (left) and tutor and pupil path length in the right, and blue lines indicating the residuals.

Point 2 - repeatability & reliability

Repeatability of path length for tutored birds

# in Alam et al's 18 bird dataset behind Fig. 2b
## load path length data
d = fread(here::here('Data/Fig_2b.csv'), header = T) %>% 
            pivot_longer( 
                cols = `1`:`15`, 
                names_to = "UMAP", 
                values_to = "path_length")  %>%   
                        data.table()
## load syllable count data
s = fread(here::here('Data/match_labels_output.csv'), header = T) 
## merge
dd = data.table(merge(d,s[,.(f2b_bird_name,superset_bird_type, superset_bird_nsyll)], all.x = TRUE, by.x = 'bird', by.y = 'f2b_bird_name')) %>% 
    setnames('superset_bird_nsyll', 'n_syll')

## repeatability for known tutored birds and with available syllable count 
m = lmer(path_length ~ (1 | bird), data = dd[superset_bird_type%in%'tutored'])
o_m = m_out(m, name = '1', dep = "path_length", save_sim = FALSE)

## repeatability controlled for syllable count
ms = lmer(path_length ~ n_syll + (1 | bird), data = dd[superset_bird_type%in%'tutored'])
o_ms = m_out(ms, name = '2', dep = "path_length", save_sim = FALSE)

# for 31 tutored birds from Fig 2c
## load data
a = fread(here::here('Data/Dat_path_length.csv'), header = T) # the dataset was created using Alam et al's Fig 2c data on 31 birds available from https://doi.org/10.18738/T8/WBQM4I/Q92O9A and using their procedure to generated 20 UMAPs - see https://github.com/ymir-k/UMAP_test/tree/main/AlamTests/3-RseedTest
a[, path_length := kdistance] #length(unique(a$bird_id)) # n = 31 birds
#ggplot(a[n_syll ==5], aes(x = bird_id, y=path_length)) + geom_point(pch = 21) + geom_line(aes(col = as.factor(iteration)))
#ggplot(a, aes(x = n_syll, y=path_length, fill = as.factor(n_syll), group=bird_id)) +geom_point(pch = 21)

## repeatability
m31 = lmer(path_length ~ (1 | bird_id), data = a)
o_m31 = m_out(model = m31, name = '3', dep = "path_length", save_sim = FALSE)

## repeatability controlled for syllable count
ms31 = lmer(path_length ~ n_syll + (1 | bird_id), data = a)
o_ms31 = m_out(model = ms31, name = '4', dep = "path_length", save_sim = FALSE)

# export
o1 = rbind(o_m,o_ms)
o1[!is.na(N), `N birds` := dd[superset_bird_type%in%'tutored', length(unique(bird))]]
o1[model %in%1, control:="none"]
o1[model %in%2, control:="for syllable count"]

o2 = rbind(o_m31,o_ms31)
o2[!is.na(N), `N birds` := length(unique(a$bird_id))]
o2[model %in%3, control:="none"]
o2[model %in%4, control:="for syllable count"]

oo = rbind(o1,o2)
oo[, `95%CI` := paste(lwr_r, upr_r, sep = "-")]
setnames(oo, c('estimate_r'), c('estimate'))
oo[response%in%'path_length', response := 'path length']
oo[effect%in%'n_syll', effect := 'syllable count']
oo[is.na(control), control :=""]
oo[effect%in%"bird_id (Intercept)", effect :="bird (Intercept)"]

oo[,.(model, control, response, N, `N birds`, type, effect, estimate,`95%CI` )]  %>%
  kbl(caption = "Table R | Repeatability of path-length in tutored birds") %>%
  kable_paper("hover", full_width = F) #%>%
Table R | Repeatability of path-length in tutored birds
model control response N N birds type effect estimate 95%CI
1 none path length 120 8 fixed (Intercept) 29.85 25.76-34.06
8 random bird (Intercept) 45% 40%-47%
8 random Residual 55% 53%-60%
2 for syllable count path length 120 8 fixed (Intercept) 15.2 6.13-23.88
8 fixed syllable count 2.74 1.17-4.35
8 random bird (Intercept) 21% 13%-28%
8 random Residual 79% 72%-87%
3 none path length 620 31 fixed (Intercept) 35.41 32.27-38.48
31 random bird (Intercept) 55% 54%-55%
31 random Residual 45% 45%-46%
4 for syllable count path length 620 31 fixed (Intercept) 9.6 2.84-16.36
31 fixed syllable count 5.33 4.01-6.68
31 random bird (Intercept) 27% 24%-30%
31 random Residual 73% 70%-76% 
  #scroll_box(width = "100%", height = "650px")

Estimates and 95% credible intervals based on 5,000 simulated values generated by the sim function from the arm R-Package. The repeatability estimates represent the effect of ‘bird’. Model 1 and 2 uses tutored birds with tracable syllable count from Alam et al’s Fig. 2b (n = 8), the models 3 and 4, their data behind Fig. 4c (n = 31)

Reliability of within pair path-length differences

# load data
w = fread(here::here('Data/Dat_path_length.csv'), header = T) # the dataset was created using Alam et al's Fig 2c data on 31 birds available from https://doi.org/10.18738/T8/WBQM4I/Q92O9A and using their procedure to generated 20 UMAPs - see https://github.com/ymir-k/UMAP_test/tree/main/AlamTests/3-RseedTest
setnames(w, c('kdistance'), c('path_length'))#w[, path_length := kdistance]
# selcet 4 syllable males
j=4
wj = w[n_syll%in%j] 
# create pair dataset (for each UMAP iteration add the path length of the second bird in a pair)
pj <- wj[, { 
  # generate unique pairs of bird IDs
  bird_pairs <- CJ(bird_a = bird_id, bird_b = bird_id, sorted = TRUE)[bird_a < bird_b]

  # merge path lengths for bird_a and bird_b
  bird_pairs <- merge(bird_pairs, .SD[, .(bird_id, path_length)], 
                      by.x = "bird_a", by.y = "bird_id", all.x = TRUE)
  setnames(bird_pairs, "path_length", "path_length_a")

  bird_pairs <- merge(bird_pairs, .SD[, .(bird_id, path_length)], 
                      by.x = "bird_b", by.y = "bird_id", all.x = TRUE)
  setnames(bird_pairs, "path_length", "path_length_b")

  bird_pairs
}, by = iteration]

pj[, path_length_diff := path_length_a-path_length_b]
pj[, delta := abs(path_length_diff)]
pj[, pr := paste(bird_a,bird_b)]

# create datasets with one iteration as a reference
bl = list() # list for the figure data
bbl = list() # list to estimate % of flipping

for(i in unique(pj$iteration)){ #pj data.table comes from the sessions above
  #i = 3
  xi = data.table(pr = pj[iteration==i, pr], pair = pj[iteration==i, round(abs(path_length_diff),4)], delta = pj[iteration==i, round(path_length_diff,4)])
  xi[pair==delta,adjust:= 1]
  xi[is.na(adjust), adjust :=-1]
  pxi = merge(pj,xi[,.(pr, pair,adjust)],all.x=TRUE)
  pxi[, delta:=path_length_diff*adjust ]# for plotting make all starting pair values from umap 5 positive
  pxi[, ref_iteration :=  i]

  # create panel titles and data to estimate % of flipping
  pxip = pxi[!iteration%in%i, if(adjust[1]==-1){sum(path_length_diff > 0)}else{sum(path_length_diff <0)}, by = list(pr,pair)] %>% setnames(old = 'V1', new = 'flip')
  pxip[, flip_per:=paste0(round(100*flip/(length(unique(pxi$iteration))-1)),'%')]  #paste0(round(100*swap/20), %)]  
  pxip[, flip_pr:=round(100*flip/(length(unique(pxi$iteration))-1))]
  pxip[, ref_iteration :=  i]

  pxi[, panel_titles := paste0("iteration: ", i, "; reclassified: ", paste0(round(mean(pxip$flip_pr)),'%'))]

  bl[[i+1]]= pxi
  bbl[[i+1]] = pxip
}

# combine (the data are then also used to estimate main text results and create ED_Fig. 2)
b = do.call(rbind,bl)
bb = do.call(rbind,bbl)

b[, panel_titles := factor(panel_titles, levels = unique(b[order(ref_iteration), panel_titles]))] # define plotting order

# estimate the flipping % mean(bb$flip_pr)

# plot
i = 3 # use the most variable iteration
pxi = b[ref_iteration%in%i] 

f2 =
ggplot(pxi, aes(x = pair, y = delta)) + 
    geom_rect(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax =0,fill = 'grey80',inherit.aes = FALSE)+
    geom_rect(xmin = -Inf, xmax = Inf, ymin = 0, ymax = Inf,fill = 'white',inherit.aes = FALSE)+
    geom_vline(aes(xintercept = pair), col = 'grey75', lwd = 0.25) +
    geom_hline(yintercept = 0, lty = 3, col = 'grey45', lwd = 0.25) +
    geom_point(alpha = 0.8, pch = 21, size = 2, fill =furt_) + 
    geom_point(data = pxi[iteration==i], aes(x = pair, y = delta), fill = orig_, pch = 21, size = 2) +
    #geom_text(data = wxp, aes(x = pair, y = -45, label = swap_per), size = 0.8/scale_size, angle = 90, hjust = 0)+
    annotate("text", x = 38.5, y = 39, label = "Initial difference", col = orig_, size = 0.9/scale_size, hjust = 1) + 
    #annotate("text", x = 39.75, y = 19, label = "Further differences", col = furt_, size = 0.9/scale_size, hjust = 1) + 
    annotate("text", x = 48, y = 0, label = "Initial long song", size = 1/scale_size, angle = 90) + 
    annotate("text", x = 49, y = 25, label = "remained the long one", size = 1/scale_size, angle = 90) + 
    annotate("text", x = 49, y = -25, label = "became the short one", size = 1/scale_size, angle = 90) + 
    labs(x = 'Initial path-length difference of a pair\n[indicated also by red dots]', y = 'Path-length difference\nfor each latent space iteration')+ #, tag = 'a'
    #scale_y_continuous(lim = c(-30.2, 30.2), breaks=seq(-30,30, by = 10))+
    scale_x_continuous(lim = c(-1,50), breaks = seq(0,50,by=10),expand = c(0,0))+
    scale_y_continuous(lim = c(-50,50), breaks = seq(-50,50,by=10), expand = c(0,0))+
    #scale_fill_manual(values = col_3) + 
    theme_bw() +
    theme(legend.position="none", 
    plot.tag = element_text(face = "bold", size = 10),
    axis.ticks = element_blank())    

ggsave(here::here("Output/rev_Fig_2_width-110mm.png"), f2, width = 17, height = 7.4*8.5/8 , unit = "cm")
f2

Figure 2 | Inconsistent long-path classification within song pairs across latent (UMAP) space iterations. The x-axis represents the initial difference in path lengths between two songs of a pair (also highlighted by red points), calculated from a single latent space iteration. The y-axis shows all within-pair path-length differences across 20 latent space iterations. Grey vertical lines highlight individual pairs. We used Alam et al.’s data on 31 tutored birds from Fig. 2c3 and generated 20 latent space iterations4. To control for variation due to syllable count, we selected 11 birds with four-syllable songs (the largest sample size for any syllable count), created 55 song pairs, and calculated within-pair path-length differences. A single latent space iteration served as the reference (“initial difference” in red), just like Alam et al. used only one iteration to define which song in a stimulus pair had the longer path. The blue points show the within-pair path-length differences for the remaining 19 iterations. The grey area highlights the 42% of cases where long-path song is reclassified as short-path in another iteration. Depicted are data for the reference iteration with the largest variation in path-length differences, but other iterations provide similar results (mean = 38% of reclasifications; Extended Data Fig. 2).

# load data
w4 = fread(here::here('Data/Dat_path_length.csv'), header = T)[n_syll%in%4] # the dataset was created using Alam et al's Fig 2c data on 31 birds available from https://doi.org/10.18738/T8/WBQM4I/Q92O9A and using their procedure to generated 20 UMAPs - see https://github.com/ymir-k/UMAP_test/tree/main/AlamTests/3-RseedTest
setnames(w4, c('kdistance'), c('path_length'))
 #summary(w[,.(bird_id, iteration, rseed, distance, path_length, n_syll)])

# estimate the true value per male (best unbiased linear predictor)
m1 = lmer(path_length ~ 1 + (1|bird_id), w4, REML=TRUE) # mixed effect model

blups_with_intercept <- data.table(coef(m1)$bird_id) %>% setnames(new = 'blup')# extracts blubs with added intercept (same as ranef(m1)$bird_id + 36.652)
blups_with_intercept[,bird_id:= unique(w4$bird_id)]

wb = merge(w4, blups_with_intercept, all.x = TRUE)

# plot the true vlaue agains the iterations
#ggplot(wb, aes(x = blup, y = path_length, fill = as.factor(blup))) + geom_point(pch = 21, col = "white", size = 4)+ labs(x = "True inherent path length",  y = "Path length for each iteration") + coord_cartesian(xlim = c(10,70), ylim = c(10,70))+ theme(legend.position="none")

# prepare the dataset for estimating the probability of incorrect assignment
samples <- data.table(t(combn(unique(wb$bird_id), 6))) # all unique subsamples of 6 birds out of 11

pp = foreach(i = unique(wb$iteration), .combine = rbind) %do% {
# i = 0
    wi = wb[iteration==i]
    
    p = foreach(j = 1:nrow(samples), .combine = rbind) %do% {
        #j = 1
        wij = wi[bird_id%in%samples[j]] 
        wij = wij[order(path_length)]
        pair = data.table(
            pair = c('#1',"#2","#3"),
            bird_a = c(wij$bird_id[2], wij$bird_id[1], wij$bird_id[3]),
            bird_b = c(wij$bird_id[5], wij$bird_id[6], wij$bird_id[4]),
            delta = c(wij$path_length[5]-wij$path_length[2],
           wij$path_length[6]- wij$path_length[1],wij$path_length[4]- wij$path_length[3]
            ),
            delta_blub = c(wij$blup[5]-wij$blup[2], wij$blup[6] - wij$blup[1],wij$blup[4] - wij$blup[3]
            ),
        set = j
        )
        return(pair)
    }
    return(p[, iteration := i])
}
#nrow(pp)

# assign 1 when the delta from one iteration is of oposite sign (long-path wrongly assigned) compare to the true difference
pp[, incorrect := ifelse((delta >= 0 & delta_blub >= 0) | (delta < 0 & delta_blub < 0), 0, 1)] 

# create a unique pair_id ensuring order is consistent
pp[, pair_ab := paste(pmin(bird_a, bird_b), pmax(bird_a, bird_b), sep = "_")]

# prepare predictions
m2 = glm(incorrect~delta, pp, family = "binomial")
bsim = sim(m2,n.sim = nsim)
v <- apply(bsim@coef, 2, quantile, prob = c(0.5))#plogis(v)
newD = data.table(delta = c(5,15,30,seq(0, max(pp$delta), length.out =100)))
newD = newD[order(delta)]
X = model.matrix(~delta, newD)
predmatrix = matrix(nrow = nrow(newD), ncol = nsim)
for(j in 1:nsim) predmatrix[,j] <-plogis(X %*% bsim@coef[j,])
newD$pred <- plogis(X %*% v) # fitted values
newD$lwr <- apply(predmatrix, 1, quantile, prob = 0.025)
newD$upr <- apply(predmatrix, 1, quantile, prob = 0.975)

# force through zero or not
if(force_0_yes==TRUE){force0=0.5-newD$pred[1]}else{force0=0}
newD$pred_0 = newD$pred + force0
newD$lwr_0 = newD$lwr + force0
newD$upr_0 = newD$upr + force0

# prepare dataset for Alam et al's pair scenario
pa =newD[delta %in% c(5, 15, 30)]
pa[, pair:=c('Pair #3', 'Pair #1', 'Pair #2')]
pa[, pred_text:=pred_0] # placement of labels
pa[pair %in% 'Pair #3', pred_text := pred_0 + 0.05]

# prepare mean values for plotting
pp[, delta_bin := cut(delta, 
    breaks = seq(min(delta), max(delta), length.out = 11), 
    include.lowest = TRUE)] #pp[, delta_bin := cut(delta, breaks = quantile(delta, probs = seq(0, 1, length.out = 11)), include.lowest = TRUE)]

summary_stats <- pp[, .(
    median_pr = median(incorrect, na.rm = TRUE), 
    mean_pr = mean(incorrect, na.rm = TRUE),
    mean_delta = mean(delta, na.rm = TRUE),
    n = length(pair)
), by = delta_bin]

# prepare densityplot
pp_0 <- pp[incorrect == 0]
pp_1 <- pp[incorrect == 1]
density_1 <- density(pp_1$delta) # Compute density for y = 1

## Convert to a dataframe
density_df_1 <- data.frame(
  x = density_1$x, 
  y = density_1$y / max(density_1$y) * 0.1 # Scale the density
)
## Flip the density so it extends *upwards* from y = 1
density_df_1$ymin <- 1 - density_df_1$y  # Bottom boundary
density_df_1$ymax <- 1                   # Top boundary (fixed at y = 1)

# prepare legend
g_leg = 
ggplot() +
  geom_point(data = summary_stats, aes(x = mean_delta, y = mean_pr, size = n), col = point_out, fill = point_fill, pch = 21, alpha = 0.8) +
  scale_size_area(
    limits = c(0, max(summary_stats$n, na.rm = TRUE)), 
    max_size = 10, 
    breaks = c(100, 1000, 5000), 
    labels = c(100, 1000, 5000), 
    guide = guide_circles(vjust = 1)
  )+
  theme_bw() +
  theme(  legend.text=element_text(size=7),
          legend.title=element_text(hjust=0.5, margin = margin(b = -1)),
          legend.key = element_rect(colour = NA, fill = NA),
          #legend.margin = margin(0,0,0,0, unit="cm"),
          legend.background = element_blank())
leg <- get_legend(g_leg)

# prepare plot
ed_f = 
ggplot() + 
  # Add marginal density plots at y = 0 and y = 1
  geom_density(data = pp_0, aes(x = delta, y = after_stat(scaled) * 0.1), alpha = 0.3, fill = "grey90", col = "grey70") + 
  geom_density(data = pp_1, aes(x = delta, y = 1 - after_stat(scaled) * 0.1),  alpha = 0.3, fill = NA, col = "grey70") +
  geom_ribbon(data = density_df_1, aes(x = x, ymin = ymin, ymax = ymax), 
              fill = "grey90", alpha = 0.3) +

  # Add means of raw data
  geom_point(data = summary_stats, aes(x = mean_delta, y = mean_pr, size = n), col = point_out, fill = "grey70", pch = 21, alpha = 0.8) +

  # prediction from the logistic regression without 95%CI, as those are pseudoreplicated        
  #geom_ribbon(data = newD,aes(ymin=lwr_0, ymax=upr_0, x=delta), fill = orig_, alpha = 0.2, show.legend = NA) +
  geom_line(data = newD,aes(x = delta, y = pred_0), col = orig_) +
  # Alam et al's pairs highlighted 
  geom_point(data = pa, aes(x = delta, y = pred_0), pch = 21, fill = "red")+
   # Add labels above the three points TODO:move Pair #3 higher
  geom_text(data = pa, aes(x = delta, y = pred_text, label = round(pred_0, 2)), vjust = -1, hjust = 0, color = "grey20", fontface = "bold", size = 1/scale_size) +  
  geom_text(data = pa, aes(x = delta, y = pred_text, label = pair), vjust = -2.5, hjust = 0, color = "grey20", size = 1/scale_size) +  
  scale_size_area(
    limits = c(0, max(summary_stats$n, na.rm = TRUE)), 
    max_size = 10, 
    breaks = c(100, 1000, 5000), 
    labels = c(100, 1000, 5000), 
    guide = guide_legend(vjust = 1)
  )+
  scale_x_continuous(expand = c(0, 0),name = "Path-length difference of a pair\n[based on one iteration]", breaks=c(0,15,30,45)) +
  scale_y_continuous(expand = c(0, 0), name = "Probability of incorrect classification\n[of long-path song within a pair]")+
  coord_cartesian(xlim =c(0,45), ylim = c(0,1))+
  #breaks = seq(0,1, by=0.25), labels = seq(0,100, by = 25)) 
  #labs(tag = "b") +
  theme_bw()+
  theme(legend.position = "none",
        axis.ticks = element_blank(),
        plot.tag = element_text(face = "bold", size = 10),
        panel.grid.major = element_blank(),  # Remove major grid lines
        panel.grid.minor = element_blank()   # Remove minor grid lines
        ) 

ed_f1 = 
ggdraw() +
  draw_plot(ggdraw(ggarrange(ed_f)) +
              draw_grob(leg, x = 0.7, y = 0.65, width = 0.3, height = 0.3),
            x = 0, y = 0, width = 1, height = 1)

ggsave(here::here("Output/rev_ED_Fig_1_width-55mm_v2.png"), ed_f1, width = 8.5, height = 8, units = "cm")  # width 8.5*.65
ed_f1

Extended Data Figure 1 | Probability of incorrect long-path classification within stimulus pairs as a function of their path-length difference in a single UMAP iteration. Density plots at the top and bottom indicate the distribution of the data. Grey points represent means for ten equally spaced x-axis intervals. The red line represents the predicted probability based on the joint posterior distribution of 5,000 simulated values from a logistic regression using the sim function from the arm R-package. The three red dots, along with their values, indicate the probability of incorrect assignment for the three stimulus pairs used by Alam et al. The data for this figure were obtained as follows. We used Alam et al.’s spectrograms on 31 tutored birds from Fig. 2c3 and generated 20 latent space iterations4. Within each iteration, for each male we calculated the shortest path length between syllable clusters. To control for variation due to syllable count, we selected 31 birds with four-syllable songs (the largest sample size for any syllable count). We fitted an intercept-only mixed-effects model with bird identity as a random intercept and using the coef() function in R extracted an unbiased estimate of the “true inherent path length” for each male, approximating an average over an infinite number of iterations. To simulate Alam et al’s playback scenario of three stimulus pair contrasting presumably “long path” against “short path” songs, for each of the 20 iterations, we created all unique subsets of six males (n = 462), each with a slightly different set of six males. Within each subset (following Alam et al’s Extended Data Fig. 5d), songs were sorted by their path length, and three song pairs created: Pair #1 contrasts the 2nd- and 5th-ranked song, Pair #2 the 1st- and 6th-ranked song, and Pair #3 the 3rd- and 4th-ranked song. This process resulted in 27720 comparisons (three pairs per subset across 462 subsets and 20 iterations). For each pair, we noted whether the song appearing as “long-path” in a single iteration was actually “short-path” based on its true inherent value. A logistic regression was fitted to estimate how the probability of incorrect classification changes with the path-length difference in a given iteration. Note that the probability that all three stimulus pairs were correctly classified is only 44% (the product of the individual correct-classification probabilities).

# uses data prepared in Figure 2 script above
# plot
ed_f2 =
ggplot(b, aes(x = pair, y = delta)) + 
    geom_rect(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax =0,fill = 'grey80',inherit.aes = FALSE)+
    geom_rect(xmin = -Inf, xmax = Inf, ymin = 0, ymax = Inf,fill = 'white',inherit.aes = FALSE)+
    geom_vline(aes(xintercept = pair), col = 'grey75', lwd = 0.25) +
    geom_hline(yintercept = 0, lty = 3, col = 'grey45', lwd = 0.25) +
    geom_point(alpha = 0.5, pch = 21, size = 2, fill =furt_) + 
    geom_point(data = b[iteration==ref_iteration], aes(x = pair, y = delta), fill = orig_, pch = 21, size = 2) +
    facet_wrap(~panel_titles, nrow = 5, ncol = 4) + 
    labs(x = 'Initial path-length difference of a pair\n[indicated also by red dots]', y = 'Path-length difference of a pair\nfor each latent space iteration')+ #subtitle = paste("Iteration #:",i))+
    #scale_y_continuous(lim = c(-30.2, 30.2), breaks=seq(-30,30, by = 10))+
    scale_x_continuous(lim = c(-1,50), ,expand = c(0,0))+ #breaks = seq(0,50,by=10)
    scale_y_continuous(lim = c(-50,50),  expand = c(0,0))+ #breaks = seq(-50,50,by=10),
    #scale_fill_manual(values = col_3) + 
    geom_text(
      data = b[panel_titles == unique(b$panel_titles)[1]][1,],  # Select first panel only
      aes(x = 40, y = -25, label = "reclassified"),
      size = 7*scale_size
    ) +
    geom_text(
      data = b[panel_titles == unique(b$panel_titles)[1]][1,],  # Select first panel only
      aes(x = 40, y = 25, label = "hold"),
      size = 7*scale_size
    ) +
    theme_bw() +
    theme(legend.position="none", 
    strip.background = element_blank(),
    axis.ticks = element_blank()
    )

# remove lables
ed_f2_g <- ggplotGrob(ed_f2) # gg$layout$name
gtable_filter_remove <- function (x, name, trim = TRUE){
  matches <- !(x$layout$name %in% name)
  x$layout <- x$layout[matches, , drop = FALSE]
  x$grobs <- x$grobs[matches]
  if (trim) 
    x <- gtable_trim(x)
  x
}
ed_f2_gg <- gtable_filter_remove(ed_f2_g, name = c("axis-b-2-5", "axis-b-4-5"), trim = FALSE) # paste0("axis-b-", c(2, 4), "-9")
ggsave(here::here(paste0("Output/rev_ED_Fig_2_width-130mm.png")), ed_f2_gg, width = 20, height = 18, unit = "cm") #width =  20*0.65
 
grid.draw(ed_f2_gg)

Extended Data Figure 2 | Inconsistent long-path classification within song pairs across latent (UMAP) space iterations. The x-axis represents the initial difference in path lengths between two songs of a pair (also highlighted by red points), calculated from a single latent space iteration. The y-axis shows all within-pair path-length differences across 20 latent space iterations. Grey vertical lines highlight individual pairs. Each panel represents a different latent space iteration serving as the reference. We used Alam et al.’s data on 31 tutored birds from Fig. 2c3 and generated 20 latent space iterations4. To control for variation due to syllable count, we selected 11 birds with four-syllable songs (the largest sample size for any syllable count), created 55 song pairs, and calculated within-pair path-length differences. A single latent space iteration served as the reference (“initial difference” in red), just like Alam et al. used only one iteration to define which song in a stimulus pair had the longer path. The blue points show the within-pair path-length differences for the remaining 19 iterations. The grey area illustrates how often the long-path song is reclassified as short-path in another iteration. The panel titles highlight the percentage of reclassifications, which was overall 38%.

# prepare data
w = fread(here::here('Data/Dat_path_length.csv'), header = T) # the dataset was created using Alam et al's Fig 2c data on 31 birds available from https://doi.org/10.18738/T8/WBQM4I/Q92O9A and using their procedure to generated 20 UMAPs - see https://github.com/ymir-k/UMAP_test/tree/main/AlamTests/3-RseedTest
setnames(w, c('kdistance'), c('path_length'))#w[, path_length := kdistance]
# selcet 4 syllable males
j=5
wj = w[n_syll%in%j] 
# create pair dataset (for each UMAP iteration add the path length of the second bird in a pair)
pj <- wj[, { 
  # generate unique pairs of bird IDs
  bird_pairs <- CJ(bird_a = bird_id, bird_b = bird_id, sorted = TRUE)[bird_a < bird_b]

  # merge path lengths for bird_a and bird_b
  bird_pairs <- merge(bird_pairs, .SD[, .(bird_id, path_length)], 
                      by.x = "bird_a", by.y = "bird_id", all.x = TRUE)
  setnames(bird_pairs, "path_length", "path_length_a")

  bird_pairs <- merge(bird_pairs, .SD[, .(bird_id, path_length)], 
                      by.x = "bird_b", by.y = "bird_id", all.x = TRUE)
  setnames(bird_pairs, "path_length", "path_length_b")

  bird_pairs
}, by = iteration]

pj[, path_length_diff := path_length_a-path_length_b]
pj[, delta := abs(path_length_diff)]
pj[, pr := paste(bird_a,bird_b)]

# create datasets with one iteration as a reference
bl5 = list() # list for the figure data
bbl5 = list() # list to estimate % of flipping

for(i in unique(pj$iteration)){ #pj data.table comes from the sessions above
  #i = 3
  xi = data.table(pr = pj[iteration==i, pr], pair = pj[iteration==i, round(abs(path_length_diff),4)], delta = pj[iteration==i, round(path_length_diff,4)])
  xi[pair==delta,adjust:= 1]
  xi[is.na(adjust), adjust :=-1]
  pxi = merge(pj,xi[,.(pr, pair,adjust)],all.x=TRUE)
  pxi[, delta:=path_length_diff*adjust ]# for plotting make all starting pair values from umap 5 positive
  pxi[, ref_iteration :=  i]

  # create panel titles and data to estimate % of flipping
  pxip = pxi[!iteration%in%i, if(adjust[1]==-1){sum(path_length_diff > 0)}else{sum(path_length_diff <0)}, by = list(pr,pair)] %>% setnames(old = 'V1', new = 'flip')
  pxip[, flip_per:=paste0(round(100*flip/(length(unique(pxi$iteration))-1)),'%')]  #paste0(round(100*swap/20), %)]  
  pxip[, flip_pr:=round(100*flip/(length(unique(pxi$iteration))-1))]
  pxip[, ref_iteration :=  i]

  pxi[, panel_titles := paste0("iteration: ", i, "; reclassified: ", paste0(round(mean(pxip$flip_pr)),'%'))]

  bl5[[i+1]]= pxi
  bbl5[[i+1]] = pxip
}

# combine (the data are then also used to estimate main text results and create ED_Fig. 2)
b5 = do.call(rbind,bl5)
bb5 = do.call(rbind,bbl5)

b5[, panel_titles := factor(panel_titles, levels = unique(b5[order(ref_iteration), panel_titles]))] # define plotting order


# plot
s_f1 =
ggplot(b5, aes(x = pair, y = delta)) + 
    geom_rect(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax =0,fill = 'grey80',inherit.aes = FALSE)+
    geom_rect(xmin = -Inf, xmax = Inf, ymin = 0, ymax = Inf,fill = 'white',inherit.aes = FALSE)+
    geom_vline(aes(xintercept = pair), col = 'grey75', lwd = 0.25) +
    geom_hline(yintercept = 0, lty = 3, col = 'grey45', lwd = 0.25) +
    geom_point(alpha = 0.5, pch = 21, size = 2, fill =furt_) + 
    geom_point(data = b5[iteration==ref_iteration], aes(x = pair, y = delta), fill = orig_, pch = 21, size = 2) +
    facet_wrap(~panel_titles, nrow = 5, ncol = 4) + 
    labs(x = 'Initial path-length difference of a pair\n[indicated also by red dots]', y = 'Path-length difference of a pair\nfor each latent space iteration')+ #subtitle = paste("Iteration #:",i))+
    #scale_y_continuous(lim = c(-30.2, 30.2), breaks=seq(-30,30, by = 10))+
    scale_x_continuous(lim = c(-1,50), ,expand = c(0,0))+ #breaks = seq(0,50,by=10)
    scale_y_continuous(lim = c(-50,50),  expand = c(0,0))+ #breaks = seq(-50,50,by=10),
    #scale_fill_manual(values = col_3) + 
    geom_text(
      data = b5[panel_titles == unique(b5$panel_titles)[1]][1,],  # Select first panel only
      aes(x = 40, y = -25, label = "reclassified"),
      size = 7*scale_size
    ) +
    geom_text(
      data = b5[panel_titles == unique(b5$panel_titles)[1]][1,],  # Select first panel only
      aes(x = 40, y = 25, label = "hold"),
      size = 7*scale_size
    ) +
    theme_bw() +
    theme(legend.position="none", 
    strip.background = element_blank(),
    axis.ticks = element_blank()
    )

# remove lables
s_f1_g <- ggplotGrob(s_f1) # gg$layout$name
gtable_filter_remove <- function (x, name, trim = TRUE){
  matches <- !(x$layout$name %in% name)
  x$layout <- x$layout[matches, , drop = FALSE]
  x$grobs <- x$grobs[matches]
  if (trim) 
    x <- gtable_trim(x)
  x
}
s_f1_gg <- gtable_filter_remove(s_f1_g, name = c("axis-b-2-5", "axis-b-4-5"), trim = FALSE) # paste0("axis-b-", c(2, 4), "-9")
ggsave(here::here(paste0("Output/rev_s_f1_gg_width-130mm.png")), s_f1_gg, width = 20, height = 18, unit = "cm") #width =  20*0.65
 
grid.draw(s_f1_gg)

Supporting Figure 1 | Inconsistent long-path classification within song pairs across latent (UMAP) space iterations using 5-syllable songs. The x-axis represents the initial difference in path lengths between two songs of a pair (also highlighted by red points), calculated from a single latent space iteration. The y-axis shows all within-pair path-length differences across 20 latent space iterations. Grey vertical lines highlight individual pairs. Each panel represents a different latent space iteration serving as the reference. We used Alam et al.’s data on 31 tutored birds from Fig. 2c3 and generated 20 latent space iterations4. To control for variation due to syllable count and mimick Alam et al.’s 5-syllable song stimuli, we selected 7 birds with five-syllable songs (the largest sample size for any syllable count), created 21 song pairs, and calculated within-pair path-length differences. A single latent space iteration served as the reference (“initial difference” in red), just like Alam et al. used only one iteration to define which song in a stimulus pair had the longer path. The blue points show the within-pair path-length differences for the remaining 19 iterations. The grey area illustrates how often the long-path song is reclassified as short-path in another iteration. The panel titles highlight the percentage of reclassifications, which was overall 41%.

Output used in the replies to the reviewer comments

gr = 
ggplot(a, aes(x= as.factor(n_syll), y = path_length)) +  
    geom_boxplot(width=0.2, position = position_nudge(-0.1), outliers = FALSE) +
     geom_dotplot(binaxis = "y", width = 0.2, position = position_nudge(x = 0.1), binwidth = 0.6, fill = point_fill, col = point_out, alpha = 0.8) + 
     labs(y ='Path length' , x ='Syllables per song [count]')+
     scale_y_continuous(lim = c(0,80), expand = c(0,0)) + 
     theme_bw()
ggsave(here::here("Output/rev_Reply_2.png"), gr, width = 8.5, height = 8.25, unit = "cm")
gr

Responses Figure 2 | Path length as a function of syllable count. Points represent song path-length for each of the 31 males and 20 iterations (n = 620). Boxplots depict median (horizontal line inside the box), the 25th and 75th percentiles (box) ± 1.5 times the interquartile range or the minimum/maximum value, whichever is smaller (bars).

Point 3 - experimental design

# load data
r = fread(here::here('Data/Fig_3c_ED_Fig_6.csv'))
r[, path_delta:=long-short]
r[, long_vs_baseline:=(trial_ed-(pre_ed+post_ed)/2)/100] # using data behind Extended Data Fig. 6, because of discrepencies in source data of Fig. 3c, see column comments for details

# baseline mean of pre and post
rx = r[, mean(long_vs_baseline), by = path_delta] %>% setnames('V1', 'mean')

f3=
ggplot() + 
    geom_hline(yintercept = 0, lty = 3, col = "grey80")+
    geom_point(aes(x = path_delta, y = long_vs_baseline, col = side_bias), data = r, cex = 2, alpha = 0.6) +
    #stat_poly_line(se = FALSE, col = col_l) +
    geom_point(aes(x = path_delta, y = mean), data = rx, pch = 23, cex = 2, fill = 'red') +
    #geom_point(aes(x = path_delta, y = mean), data = rx, pch = 21, cex = 1, fill = 'red', alpha = 0.6) +
    
    annotate("text", x = 17, y = 0.26, label = "Means", size =2.75, hjust = 0, col = 'red') + 
    annotate("text", x = 10, y = 0.55, label = "Long-path song played in:", size =2.75, hjust = 0) + 
    annotate("text", x = 13, y = 0.515, label = "preferred arm", size =2.75, col = col_[2], hjust = 0) + 
    annotate("text", x = 13, y = 0.48, label = "non-preferred arm", size =2.75, col = col_[1], hjust = 0) + 

   # stat_cor(cor.coef.name = "r", aes(label = after_stat(r.label)),  col = col_R, r.accuracy = 0.1, label.x = 17.5, label.y = 0.6-(0.6+.03)*per_*1.3,hjust = 0.5, vjust = 0, cex = 3) +

    coord_cartesian(xlim = c(0, 35), ylim = c(-.09, .6)) +
    scale_x_continuous(breaks = seq(0,30, by=10), expand = c(0,0))+
    scale_y_continuous(breaks = seq(0,.6, by=.2), expand = c(0,0))+
    scale_color_manual(guide = "none", values = col_) +
    
    annotate("text", x = 15, y = -0.06, label = "Pair #1", size =3) + 
    annotate("text", x = 30.29, y = -0.06, label = "Pair #2", size =3) + 
    annotate("text", x = 5.72, y = -0.06, label = "Pair #3", size =3) + 
    #scale_shape_manual(guide = "none", values = c(21, 22, 23)) +
    #guides(fill = guide_legend(override.aes=list(shape=21))) +
    labs(x = 'Path-length difference of a pair\n[long-path minus short-path song]\n', y = 'Change in proportion of time\nwith long-path song\n[trial minus baseline]')+
    theme_bw()+
    theme(axis.ticks = element_blank(),
        panel.grid.major = element_blank(),  # Remove major grid lines
        panel.grid.minor = element_blank())#,text = element_text(family = "Arial Unicode MS"))

ggsave(here::here("Output/rev_Fig_3_width-55mm_v3.png"), f3, width = 8.5, height = 8.3, unit = "cm")#7.4*8.5/8, 

f3

Figure 3 | Female preference for the presumed long-path song in relation to path-length differences within stimulus pairs. The x-axis represents the difference in path lengths between the two songs of each stimulus pair. Notably, only stimulus pair #2 (at the right end of the plot) shows a strong contrast between a very short and a very long path, yet it elicited the weakest female preference for the long-path song. The y-axis quantifies the strength of female preference for the long-path song, calculated as the change in the proportion of time spent in the choice-chamber arm with the long-path playback during the trial compared to baseline (the mean of pre- and post-trial periods). Dots represent individual observations for each stimulus pair, with dot colour indicating whether the long-path song was played in the arm preferred by the female during the pre-trial period (orange) or in the other arm (blue), illustrating the unbalanced arm-assignment. Red diamonds represent the mean female response to each stimulus pair, indicating the effective sample size of three.

Session information

The figures were generated using the below indicated version of R and its related packages, running on the below indicated macOS. The installation of R and related R-packages takes only a few mintues. How to install R and R-packages is described here and here

pander(sessionInfo()) # pander(sessionInfo(), compact = FALSE)

R version 4.2.0 (2022-04-22)

Platform: x86_64-apple-darwin17.0 (64-bit)

locale: en_GB.UTF-8||en_GB.UTF-8||en_GB.UTF-8||C||en_GB.UTF-8||en_GB.UTF-8

attached base packages: grid, stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: tidyr(v.1.3.1), rptR(v.0.9.22), patchwork(v.1.3.0), pander(v.0.6.5), legendry(v.0.2.0), kableExtra(v.1.4.0), gridExtra(v.2.3), grafify(v.4.0.1), ggpubr(v.0.6.0), ggpmisc(v.0.5.4-1), ggpp(v.0.5.5), ggplot2(v.3.5.1), foreach(v.1.5.2), data.table(v.1.16.2), cowplot(v.1.1.3), arm(v.1.13-1), lme4(v.1.1-35.5), Matrix(v.1.6-2), MASS(v.7.3-60) and rmarkdown(v.2.29)

loaded via a namespace (and not attached): nlme(v.3.1-163), rprojroot(v.2.0.4), numDeriv(v.2016.8-1.1), tools(v.4.2.0), backports(v.1.5.0), bslib(v.0.8.0), utf8(v.1.2.4), R6(v.2.5.1), rpart(v.4.1.21), Hmisc(v.5.2-0), mgcv(v.1.9-0), colorspace(v.2.1-1), nnet(v.7.3-19), withr(v.3.0.2), tidyselect(v.1.2.1), emmeans(v.1.10.7), compiler(v.4.2.0), textshaping(v.0.3.7), cli(v.3.6.3), quantreg(v.5.99.1), htmlTable(v.2.4.3), SparseM(v.1.84-2), xml2(v.1.3.6), sandwich(v.3.0-2), labeling(v.0.4.3), sass(v.0.4.9), scales(v.1.3.0), checkmate(v.2.3.2), mvtnorm(v.1.3-2), systemfonts(v.1.0.5), stringr(v.1.5.1), digest(v.0.6.37), foreign(v.0.8-85), minqa(v.1.2.8), svglite(v.2.1.2), base64enc(v.0.1-3), pkgconfig(v.2.0.3), htmltools(v.0.5.8.1), fastmap(v.1.2.0), htmlwidgets(v.1.6.4), rlang(v.1.1.4), rstudioapi(v.0.17.1), jquerylib(v.0.1.4), generics(v.0.1.3), farver(v.2.1.2), zoo(v.1.8-12), jsonlite(v.1.8.9), dplyr(v.1.1.4), car(v.3.1-3), magrittr(v.2.0.3), polynom(v.1.4-1), Formula(v.1.2-5), Rcpp(v.1.0.13-1), munsell(v.0.5.1), fansi(v.1.0.6), abind(v.1.4-8), lifecycle(v.1.0.4), stringi(v.1.8.4), multcomp(v.1.4-25), yaml(v.2.3.10), carData(v.3.0-5), lattice(v.0.22-5), splines(v.4.2.0), knitr(v.1.49), pillar(v.1.9.0), boot(v.1.3-28.1), estimability(v.1.5.1), ggsignif(v.0.6.4), codetools(v.0.2-19), glue(v.1.8.0), evaluate(v.1.0.1), vctrs(v.0.6.5), nloptr(v.2.1.1), MatrixModels(v.0.5-3), gtable(v.0.3.6), purrr(v.1.0.2), cachem(v.1.1.0), xfun(v.0.49), xtable(v.1.8-4), broom(v.1.0.7), coda(v.0.19-4), rstatix(v.0.7.2), ragg(v.1.2.6), viridisLite(v.0.4.2), survival(v.3.5-7), tibble(v.3.2.1), lmerTest(v.3.1-3), iterators(v.1.0.14), cluster(v.2.1.4), TH.data(v.1.1-2) and here(v.1.0.1)

Demo & Reproduction Instructions

We do not provide any demo data, because the three csv files used in the analyses are small and already contain all the data used in the analyses. To generate the outputs, download the files into Data folder of your project’s root directory and create a folder Output. Open R software and the Fig_point_1.r or Fig_point_2.r R-script (for outputs on initial submission) and SI.R for the final outputs, install the R-packages indicated under the tools by running install.packages(c('arm','cowplot', 'data.table', 'foreach','ggplot2', 'ggpmisc', 'ggpubr', 'grafify', 'grid', 'gridExtra, 'kableExtra', 'legendry', 'pander', 'patchwork','rptR','tidyr') and run the scripts. The pngs of the figures and the html document of the whole supporting information will be generated in Output.